! module *polynomial_basis*
!
! Purpose:
!   (1) manage polynomial basis
!   (2) provide polynomial information
!
! Author:
!   Name: DONG Li
!   Email: dongli@lasg.iap.ac.cn

module polynomial_basis

    use message
    use run_manager

    implicit none

    private

    public polynomial_basis_init
    public polynomial_basis_evaluate
    public polynomial_basis_derivative
    public polynomials

    type polynomial
        integer :: num_monomial
        integer, allocatable :: exponents(:,:)
    end type polynomial

    integer :: num_dim

    type(polynomial), allocatable :: polynomials(:)

contains

    subroutine polynomial_basis_init(num_dim_, max_exponent)
        integer, intent(in) :: num_dim_, max_exponent

        integer :: i, j ,k, l, m, n
        integer :: num_row
        integer, allocatable :: ranges(:,:)

        num_dim = num_dim_

        allocate(polynomials(max_exponent))
        do l = 1, max_exponent
            num_row = l+1
            allocate(ranges(2:num_row,num_dim))
            ranges(2,:) = 1
            do i = 3, num_row
                do j = 1, num_dim
                    ranges(i,j) = sum(ranges(i-1,j:num_dim))
                end do
            end do
            ! total number of monomials
            polynomials(l)%num_monomial = sum(ranges)+1
            allocate(polynomials(l)%exponents&
                (polynomials(l)%num_monomial,num_dim))
            polynomials(l)%exponents(1,:) = 0
            do i = 2, num_dim+1
                polynomials(l)%exponents(i,:) = 0
                polynomials(l)%exponents(i,i-1) = 1
            end do
            m = num_dim+2 ! offset for exponents
            do i = 3, num_row
                n = 1 ! offset base for exponents in previous row
                do k = 2, i-2
                    n = n+sum(ranges(k,:))
                end do
                do j = 1, num_dim ! loop for ranges in current row
                    do k = 1, sum(ranges(i-1,j:num_dim)) ! loop for monomial
                        polynomials(l)%exponents(m,:) = &
                            polynomials(l)%exponents(n+k,:)
                        polynomials(l)%exponents(m,j) = &
                            polynomials(l)%exponents(m,j)+1
                        m = m+1
                    end do
                    n = n+sum(ranges(i-1,j:j))
                end do
            end do
            deallocate(ranges)
        end do

        return
    end subroutine polynomial_basis_init

    real function polynomial_basis_evaluate(x, order, c) result(res)
        real, intent(in) :: x(num_dim)
        integer, intent(in) :: order
        real, intent(in) :: c(*)

        integer :: i, j
        real :: temp

        res = 0.0
        do j = 1, polynomials(order)%num_monomial
            temp = 1.0
            do i = 1, num_dim
                temp = temp*x(i)**polynomials(order)%exponents(j,i)
            end do
            res = res+c(j)*temp
        end do
    
        return
    end function polynomial_basis_evaluate

    real function polynomial_basis_derivative(x, order, c, l, o) result(d)
        real, intent(in) :: x(num_dim)
        integer, intent(in) :: order
        real, intent(in) :: c(*)
        integer, intent(in) :: l
        integer, intent(in) :: o

        integer :: i, j
        real :: temp

        if(o > order) then
            call message_show(ERROR_MSG,       &
                "polynomial_basis",            &
                "polynomial_basis_derivative", &
                "Differentiation order is "    &
                // "larger than maximum "      &
                // "exponent of polynomial")
            call run_manager_end_run
        end if

        d = 0.0
        do j = 1, polynomials(order)%num_monomial
            temp = 1.0
            do i = 1, num_dim
                if(i == l) then
                    if(polynomials(order)%exponents(j,i) /= 0) then
                        temp = temp*polynomials(order)%exponents(j,i) &
                            *x(i)**(polynomials(order)%exponents(j,i)-1)
                    else
                        temp = 0.0
                        exit
                    end if
                else
                    temp = temp*x(i)**polynomials(order)%exponents(j,i)
                end if
            end do
            d = d+c(j)*temp
        end do

        return
    end function polynomial_basis_derivative

end module polynomial_basis
